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Abstract 


In  considering  potential  replacements  for  computer  facilities,  such  as  that  for  the  DREA 
DECSYSTEM  20/60,  comparison  standards  must  be  specified.  In  the  comparison  standards, 
two  notable  characteristics  that  should  be  specified  are  the  speed  and  accuracy  of  numerical 
computations.  HILBRT  is  a  benchmarking  program  which  may  be  used  for  this  purpose.  It 
is  based  on  the  inversion  of  finite  segments  of  the  infinite  Hilbert  matrix.  The  theory,  program, 
and  example  results  are  presented.  < 


R6sum6 


L'examen  des  systEmes  de  replacement  potentiels  pour  les  systEmes 
informatiques  comme,  par  example,  le  DECSYSTEM  20/60  du  CRDA,  exige 
l’Etablissement  de  critEres  de  comparaison.  II  est  nEcessaire  de  prEciser 
deux  caractEristiques  importantes  dans  ces  critEres  de  comparaison, 
c'est-E-dire  la  vitesse  et  1* exactitude  des  calculs  numEriques.  Le  programme 
d’ Evaluation  des  performances  HLBERT  peut  servir  E  cette  fin.  II  est  fondE 
sur  1* inversion  de  segments  finis  de  la  matrice  de  Hilbert  infinie.  Le  texte 
expose  les  principes  thEoriques,  dEcrit  le  logiciel  et  donne  un  exemple  des 
rEsultats  obtenus. 
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1  Introduction 


In  order  to  make  an  objective  comparison  between  computers,  one  must  consider  many 
features.  One  may  consider  such  things  as  user  friendliness,  software  availability,  hardware 
availability,  speed,  and  accuracy.  This  technical  communication  describes  one  method  of 
measuring  and  comparing  the  features  of  computer  speed  and  accuracy.  The  method  is  based 
on  the  inversion  of  finite  segments  of  the  Hilbert  matrix. 

The  Hilbert  matrix  is  a  classic  matrix  in  numerical  analysis.  As  the  finite  segment  size 
increases,  the  segments  become  increasingly  ill-conditioned.  HILBRT  is  a  benchmark  program 
written  in  ANSI  Fortran  ‘77  which  inverts  finite  segments  ot  the  infinite  Hilbert  matrix  using 
IMSL  Version  9.2[IMSL,  1984],  then  computes  accuracy  measurements  of  the  solution.  Two 
measures  are  used.  The  first  is  to  compute  the  sup-norm  of  the  difference  between  the  identity 
matrix  and  the  product  of  the  Hilbert  matrix  and  its  IMSL  approximated  inverse.  The  second 
measure  is  the  sup-norm  of  the  difference  between  the  IMSL  computed  inverse  and  a  computed 
closed  form  solution.  The  results  are  then  written  on  the  file  “HIL.DAT.”  Output  consists 
of  the  matrix  order  and  both  sup-norms  for  each  matrix.  These  are  calculated  for  matrices 
starting  with  order  one  and  increasing  by  one  each  time.  This  is  repeated  until  the  IMSL 
subroutine  is  unable  to  find  the  solution  or  an  order  of  thirty  is  reached.  Calculations  for  those 
matrices  of  orders  one  through  four  are  repeated  12000  times.  The  CPU  time  required  for  the 
repetitions  is  calculated  and  written  on  the  output  file  to  allow  comparison  of  systems.  The 
total  CPU  time  required  is  also  included. 

In  the  technical  specifications  for  the  replacement  of  DREA’s  DECSYSTEM  20/60,  the 
IMSL  library  is  included  as  an  essential  requirement  so  that  the  use  of  this  library  in  the 
benchmark  program  HILBRT  is  justified.  In  essence,  the  performance  characteristics  of  speed 
and  accuracy  depend  on  the  performance  of  the  subroutines  in  the  IMSL  library;  and,  in  turn, 
the  performance  of  the  subroutines  depends  on  the  computer  wordlength  and  the  floating-point 
format  being  used.  As  a  result,  this  benchmark  measures  speed  and  accuracy  in  an  IMSL 
environment. 


2  Operational  Section 


2.1  Identification 

PROGRAM  HILBRT 

A  Fortran  program  to  measure  computer  speed  and  accuracy.  It  will  be  one  of  many 
benchmark  programs  used  to  assist  in  the  selection  of  a  new  DREA  central  computing  facility. 

Author:  James  A.  Theriault 

2.2  Hardware/Software  Environment 

HILBRT  requires  the  following  software: 

-  ANSI  Fortran  ‘77 

•  Fortran  subroutine  LINV2F  from  IMSL  Version  9.2 

-  Fortran  subroutine  LEQT2F  from  IMSL  Version  9.2 

-  Fortran  subroutine  CPUTIM  from  DEC  20/60  library.  This  may  be  replaced  by  similar 
routines  on  other  computers. 

HILBRT  may  be  run  in  batch  mode.  There  is  no  user  input.  All  output  is  directed  to 
a  named  disk  file,  “HIL.DAT.” 

2.3  Purpose 

The  purpose  of  HILBRT  is  to  benchmark  computers  in  an  IMSL  environment  for  the  DEC 
20/60  replacement.  It  does  so  by  inverting  finite  segments  of  the  infinite  Hilbert  matrix. 


where 


(2.2) 


(n)  1 

a-  /  = - . 

,J  i  +  j-l 

Its  inverse  may  be  given  in  a  closed  form  as  shown  in  section  3.1.  The  inverse  may  also  be 
found  using  various  numerical  techniques  such  as  the  Crout  algorithm  for  Gaussian  elimination 
utilized  by  the  IMSL  routine,  LEQT2F. 

HILBRT  uses  both  methods  to  compute  the  inverse,  then  compares  the  results  for  accuracy. 


2.5  Input 

The  program  requires  no  user  input. 


2.6  Output 

The  following  is  the  output  file,  ‘HIL.DAT,’  from  the  DREA  DEC  20/60. 


ORDER 

NORM  OF  A+D-I 

NORM  OF  D-] 

1 

. OOOOOE+OO 

.OOOOOE+OO 

2 

. OOOOOE+OO 

. 23842E-06 

3 

.47684E-06 

. 85831E-04 

4 

. 76294E-05 

. 13245E-01 

6 

•24414E-03 

51.516 

6 

.78125E-02 

31173. 

7 

. 18750 

. 17929E+08 

MATRIX  SIZE  8  UNABLE  TO  BE  INVERTED.  IMSL  ERROR  *  120 

TIME  AFTER  12000  REPETITIONS  OF  THE  FIRST  FOUR  ORDERS  -  4.984 
TOTAL  CPU  TIME  -  4.991  MINUTES 


A  refers  to  the  finite  segment  of  the  Hilbert  matrix.  I  refers  to  the  identity  matrix.  D 
and  B  refer  to  the  IMSL  calculated  approximation  and  the  closed  form  approximation  of  the 
inverse  respectively. 


2.7  Operation 

The  sequence  of  monitor  instructions  for  executing  HILBRT  on  the  DREA  DEC  20/60  is 
given  by 


•COMPILE  DREA : <MAG . THERIAULT>HILBRT . FOR 

•LOAD  DREA : <MAG . THERI AULT>HILBRT ,  <LIB>DECLIB/LIB ,  <LIB>IMSL/LIB 
•SAVE 

CRUM  HILBRT 

A  similar  sequence  of  commands  should  be  available  on  all  benchmarked  computers.  Care 
should  be  taken  to  ensure  that  the  calculated  CPU  time  does  not  include  time  for  compile  and 
load  operations.  Compiler  optimization  must  not  be  used. 

HILBRT  requires  no  further  user  input. 

2.8  Restrictions  &  Limitations 

A  maximum  size  for  the  finite  segments  of  the  Hilbert  matrix  is  30  by  30.  This  may  be 
changed  by  modifying  the  value  of  NMAX  in  the  parameter  statement  of  the  program.  It  is 
assumed  that  the  computer  can  handle  Hilbert  matrix  segments  up  to  order  n  =  4. 

2.9  Subroutines  Called 

The  following  subroutines  and  functions  are  called  by  HILBRT. 

CMPINV  -  Subroutine  to  compute  inverse  of  Hilbert  matrix  and  appropriate  norms.  It 
has  been  supplied  with  the  program. 

PERM  -  Real  function  which  computes  a  permutation.  It  has  been  supplied  with  the 
program. 

CPUTIM  -  A  DEC  20/60  dependent  subroutine  to  return  the  CPU  time  in  milliseconds 
since  the  start  of  execution.  Systems  to  be  compared  should  have  a  similar  subroutine. 

LINV2F  -  IMSL  version  9.2  subroutine  to  invert  matrix.  Calls  IMSL  subroutine  LEQT2F 
to  perform  inversion. 

2.10  Memory  Requirements 

The  following  is  a  list  of  the  files  involved  and  their  sizes  on  the  DREA  DEC  20.  HIL.DAT 
is  the  output  file.  HILBRT. FOR  is  the  Fortran  source  for  HILBRT.  HILBRT.REL  and 
HILBRT.EXE  are  the  relocatable  and  the  saved  executable  versions  of  HILBRT.  IMSL.REL 
is  the  IMSL  library. 

FILE  Characters(7  bit) 

DREA :<MAG.THERIAULT>HIL. DAT  457 

DREA : <MAG . THERIAULT>HILBRT . FOR  2958 
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yw 


FILE  Words (36  bit) 

DREA : <MAG . THERIAULT>HILBRT . EXE  8704 

DREA : <MAG . THERIAULT>HILBRT . REL  1427 

PS : <LIB>IMSL . REL  400267 


3  Technical  Section 


3.1  Details  of  Working 


The  finite  segment  of  order  n  of  the  infinite  Hilbert  matrix  is  given  by  [Gregory  et  al,  1969]: 


An  =  [-l”’! 

(3.1) 

where 

«M  =  1 

i+j-1 

(3.2) 

for  1  <  «  <  n  and  1  <  j  <  n. 

Its  inverse  may  be  found  numerically.  It  is  also  available  in 

the  following  closed  form 

[Gregory  et  al,  1969]. 

a* 

1! 

* 

(3.3) 

=  i»!?i, 

(3.4) 

where 

i(n)  _  #_! y+j _ (w  +  i-  l)l(n  +  j-  1)! _ 

*'J'  1  }  {i  +  3  ~  1)  ((<■ -1)!(J -l)0’(n  -•■)!(«■ -  j)\ 


(3.5) 


This  can  be  computed  for  successive  values  of  n  by  using  the  recursion, 


fc(n+i)  =  i(n)  (  n  +  i  \  /  n  +  j  \ 
m  «\n-i  +  l/U-J  +  l/ 

for  1  <  i  <  n  and  1  <  j  <  n.  Additional  elements,  may  be  computed  as  follows: 


°n+lj  ~  °j,n+l 


=  /  iVhW  (2n  +  l)!/n!  W(n  +  j  -  l)!/n!\ 

K  ’  V(j-l)!(n-j  +  l)!A  0  -  1)1  J 


(3.6) 


(37) 

(38) 


for  1  <  j  <  n  +  1. 

HILBRT  uses  this  and  the  IMSL  subroutine  LEQT2F  to  compute  the  inverse.  The  two 
methods  are  compared  and  the  results  of  the  comparison  are  returned  as  output. 

Define  the  sup-norm,  Halloo  =  sup,4  |c«jl  =  ■axlj|c,j|.  Let  An  be  defined  to  be  a  finite 
segment  of  the  Hilbert  matrix  as  above.  Let  Dn  be  the  computed  approximation  of  the  inverse 
using  LEQT2F  and  Bn  be  the  inverse  computed  using  the  above  closed  form. 


The  first  error  measure  is  the  sup-aorm  of  the  difference  between  the  identity  matrix  and 
the  product  of  the  Hilbert  matrix  An  and  its  IMSL  approximated  inverse  Dn: 


€{in)  =  ||  AnDn-In 

=  Halloo 


=  max  |ejj| 


(3.9) 

(3.10) 

(3.11) 


where  /„  is  the  identity  matrix  of  order  n. 

The  second  error  is  the  sup-norm  of  the  difference  between  the  IMSL  approximated  inverse 
Dn  and  the  closed  form  inverse  Bn- 


=  !Pn-£n||c 
=  nice 

=  max  |u,j|. 


(3.12) 

(3.13) 

(3.14) 


The  computation  of  the  inverses  and  comparisons  are  repeated  12000  times  for  matrices  of 
orders  one  through  four.  After  this  the  computations  continue  for  increasing  n  until  the  IMSL 
subroutine  is  unable  to  compute  the  inverse  or  until  a  maximum  order  of  thirty  is  reached. 

3.2  Flowchart 

Figure  3.1  show  a  flowchart  of  the  program  HILBRT. 

A  DEC  20/60  dependent  subroutine,  CPUTIM,  is  utilized  to  return  the  CPU  time  in 
milliseconds.  CPUTIM  is  called  at  the  start  of  HILBRT  to  return  the  CPU  time  since  the 
start  of  the  job.  It  is  called  again  after  the  first  four  matrices  have  been  inverted  12000  times 
to  return  the  time  after  execution.  It  is  called  a  final  time  after  all  calculations  are  complete 
in  order  to  compute  the  total  CPU  time.  The  elapsed  times  are  then  calculated  in  minutes 
and  written  on  the  output  file. 

A  subroutine  similar  to  CPUTIM  is  usually  available  on  most  computers. 


3.3  Important  Variables 


M 

NMAX 

A(NMAX.NMAX) 
B( NMAX, NMAX) 
D( NMAX, NMAX) 
N0RMI 


matrix  order 

matrix  order  (=  N) 

maximum  matrix  order 

finite  segment  of  Hilbert  matrix 

inverse  of  A  computed  from  closed  form  solution 

inverse  of  A  computed  by  LEQT2F 

II  A  D  -  I  IU 


NORM}  ||  D  -  B  | |oo 

HIM  time  at  start  of  job  (in  msec) 

ITIM1  time  at  end  of  job  (in  msec) 

ITIM2  time  after  12000  repetitions  of  first  4  matrices  (in  msec) 


4  Summary 


HILBRT  is  a  program  for  benchmarking  computers.  It  has  been  written  in  ANSI  Fortran 
‘77  with  the  addition  of  a  call  to  a  DEC  20/60  library  subroutine  to  determine  the  CPU  time.  It 
also  requires  a  number  of  IMSL  Version  9.2  subroutines. 

The  program  creates  and  inverts  finite  segments  of  the  infinite  Hilbert  matrix  of  orders 
one  through  to  the  largest  possible  on  the  particular  computer.  The  program  terminates  when 
an  IMSL  error  is  detected  or  when  order  30  is  completed.  The  errors  in  approximating  the 
inverses  are  written  on  the  disk  file  “HIL.DAT"  along  with  the  time  required. 
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Appendix  A  Program  Listing 


PROGRAM  HILBRT 

INVERT  FINITE  SEGMENT  OF  INFINITE  HILBERT  MATRIX 

LAST  REVISION:  DECEMBER  1.  1080 
AUTHOR:  J.  A.  THERIAULT 

PARAMETER  NMAX-30.IRPT  -  12000 

REAL  B(NMAX.NMAX) ,A(NMAX ,NMAX) .D(NMAX.NMAX) .NORMI .NORMD 

CALL  TO  DEC20  DEPENDENT  SUBROUTINE  (ITIM  -  TIME  FROM  START  IN  MSEC) 

CALL  CPUTIM(ITIM) 

OPEN (UNIT  -  22.  FILE  «  'HIL.DAT') 

WRITE (22 ,100) 

DO  8  II  -  l.IRPT 
B(l.l)  -  1. 

DO  8  N  -  1.4 
M  -  N 

CALL  CMPINV(A,B.D,M.NMAX, NORMI, NORMD. IER) 

IF  (II.EQ.l)  WRITE (22, 102)  N.NORMI, NORMD 
IF  (IER.NE.O)  GOTO  0 
CONTINUE 

IF  (IER.NE.O)  THEN 
WRITE(22.105)  I 
STOP 
ENDIF 

CALL  TO  DEC20  DEPENDENT  SUBROUTINE  (ITIM2  -  TOTAL  ELAPSED  TIME  IN  MSEC) 
CALL  CPUTIM(ITIM2) 

N  -  4 
N  -  N+l 
M  -  N 

CALL  CMP INV(A.B.D,M.NMAX, NORMI, NORMD .IER) 

IF  (IER.NE.O)  GOTO  200 
WRITE(22,102)  N.NORMI.  NORMD 
IF  (N.LE.NMAX)  GOTO  11 


tots 

1  * 

,  ■*,-  '  ■  "  ,  »*  •.'<  V  V.  .*  ;■  •  i-  | 

.•  V 


,IVI] 

<■'•3 


m 

•w 

Mi 

■  O'* 


m 


i 

m 

m 

* 
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WRITE (*.104)  N.IER 
WRITE (22. 104)  N.IER 

CALL  TO  DEC20  DEPENDENT  SUBROUTINE  (ITIM1  -  TOTAL  ELAPSED  TIME  IN  MSEC) 
CALL  CPUTIM(ITIMl) 

WRITE (22. 100)  IRPT, (ITIM2-ITIM)/60000. 

WRITE (22. 103)  (ITIMl-ITIM)/00000. 

FORMAT (/  ’  ORDER  NORM  OF  A*D-I  NORM  OF  D-B’) 

FORMAT  (I4.0X.610.6.4X.G10.B) 

FORMAT ( ’  TOTAL  CPU  TIME  -#.F.  *  MINUTES') 

FORMAT  (/’  MATRIX  SIZE  ‘.14.*  UNABLE  TO  BE  INVERTED.’. 

/’  IMSL  ERROR  14./) 

FORMAT (’  ERROR  WITH  ORDER  LESS  THAN  4’.  I) 

FORMAT ( *  TIME  AFTER  ’.17.'  REPETITIONS  OF  THE  FIRST’, 

/’  FOUR  ORDERS  .F) 

STOP 

END 

SUBROUTINE  CMPINV ( A . B . D . N , NMAX . NORMI . NORMD . IER) 

REAL  B(NMAX, NMAX) .A(NMAX. NMAX) .D(NNAX. NMAX) .NORMI. NORMD 
REAL  WORK (1000) 

NORMI  -  0. 

NORMD  -  0. 

DO  2  I  -  l.N-1 

A(I.N)  -  l./U+V-l.) 

A(N.I)  -  A(I.N) 

A(N.N)  -  l./(2*N-l.) 

IDGT  -  0 

CALL  LINV2F ( A. N. NMAX, D, IDGT, WORK, IER) 

IF  (IER.NE.O)  GOTO  8 
DO  4  I  -  l.N 
DO  4  J  -  l.N 
C  •  0. 

DO  3  K  -  l.N 

C  -  C  ♦  A(I,K)*D(K, J) 

IF  ((I.NE. J) .AND. (ABS(C) .GT. NORMI))  NORMI  -  ABS(C) 

IF  ((I.EQ.J).AND.(ABS(C-1.).GT.N0RMI))  N0RMI-ABS(C-1.) 

IF  (ABS(D(I,J)-B(I,J)).GT. NORMD)  NORMD-ABS (D ( I . J)-B(I, J)) 
CONTINUE 

DO  5  I  -  l.N 
DO  5  J-l.N 


B(I.J)  -  B(I, J)*((N*J)/(N+1 . -J))*((N+I)/(M+1.-I)) 

DO  0  J  -  l.N+1 

TMP  •  PEHM(J-l.l) 

B(N+1 , J)  -  (-1)**(N*J-1)*(PERM(2*N*1 ,H)/TMP/PERM(N*1- J.l)) 
*  *(PERM(N+J-1,N)/TMP) 

B(J.N-H)  -  B(N+1 , J) 

CONTINUE 

RETURN 

END 

REAL  FUNCTION  PERM(M.N) 

COMPUTE  MI/NI  WHERE  M  >  N 
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